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All forecasts are wrony. As a result it is rarely adequate to .jive a 
forecast by itself. Some measure of the reliability or accuracy of the fore- 
cast is also required. Such a measure is usually provided by a confidence 
interval. This is an interval within which the future value we are torecast- 
iny will fall with a prescribed probability. The philosophy and construction 
of such intervals are well documented in a number of forecasting texts and we 
will not consider them further here. See, for example, Brown (1962), 

Montgomery and Johnson (1976) and Bowerman and O'Connell (1979). 

We simply note that an essential ingredient in all the expressions and 
formulae involved in such intervals is the variance of the forecast error. 
Consider the T-step ahead forecast made at time t , i.e. the forecast of 

^ /S 

made at time t , X^(T) say. The corresponding error is e^(T) = 

2 2 

We shall denote its variance, var[e^(T)], by a-j- . It is clear that Oj must 
play a central role in any discussion of the usefulness of X^(T) in predicting 

^t+T * 

2 

The value of o-j- can usually be computed on the assumption that the under- 
lying model is valid. By underlying model we mean the model which the forecast- 
ing system assumes is generating the data. The variance can certainly be derived 
for the ARMA models of Box and Jenkins, and the interested reader is directed to 
their book (1976). It can also be derived for all the Exponential Smoothing 
systems discussed by Brown (1962) and others considered more recently by 
Sweet (1981). 

However, no results are available for Winters' seasonal systems. This is 
particularly unfortunate since these systems are amongst the most commonly used 
in practice. They are relatively simple to implement and intuitively appealing. 
Both systems are described in detail in most forecasting texts, e.g. Montgomery 
and Johnson (1976), Thomopoulos (1980), and Bowerman and O'Connell (1979). In 
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this last work, approximate confidence intervals are given for Viinters' systems. 
As we shall see, however, these intervals poorly reflect the behaviour of the 
error variance. 

The Winters' Multiplicative seasonal form is a non-linear system and it is 
difficult to see how any useful information about the forecast error can be 
obtained directly. Approximation and simulation appear to be the most sensible 
tools here. In the case of the additive model, however, we can obtain some 
results . 

The purpose of this paper is to give a procedure for deriving the variance 
of the forecast error for the Winters' Additive Seasonal forecasting system. We 
deal with not only point forecasts of future values, > l^ut also cumulative 

forecasts, i.e. forecasts of Y , = ^ X . , made at time t . 

Expressions are given for the variances of the individual components of the 
forecasts, i.e. level, trend gradient and seasonal factors, and for the covari- 
ances between them. This allows the construction of the variance of any par- 
ticular forecast and hence of the corresponding forecast error. It also enables 
us to construct confidence intervals for the components themselves, and linear 
combinations of them. 

Two special cases are considered, viz. when there is no trend, and when 
there is no seasonality. This latter is the well-known Holt-Winters non-seasonal 
forecasting system. In both cases, closed expressions are obtained for the 
variances. In the completely general seasonal case, the expressions derived 
involve a number of unknowns which are obtained by solving a set of linear 
equations. 

Some necessary discussion is also given about the stability of these 
forecasting systems. 
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1. THE RESULTS 



The underlying model assumed by the Winters' additive seasonal system is 
given by 



= m + bt + S|^ + a^ (1) 

where t = rn + k , and n is the length of the season. The values 

{S|^ : k = 1,2,. ..,n} are the additive seasonal factors and {a^} is a sequence 

of independent identically distributed random variates of zero mean and 
2 

variance o . 

The standard form of the corresponding forecasting system is given by 



"’t “O^^t'^t-n^ ^t-1^ 



The T-step ahead forecast is 

Xt(T) = m^ + Tb^ + 

where T = rn + k, (k = l,2,...,n; r ^ 0) . 

The cumulative forecast, i.e. the forecast of 



( 2 ) 



(3) 



T . T . 

''t.T ' Tj *t+i ''t,T = i, • 

Our derivation of the forecast error variances parallels the development 

given for Exponential Smoothing models by Brown (1962). By definition, X^{T) 

depends on only past values of a^ and so is independent of • Thus, 

al = + var[X.{T)J . Further, X (T) can be expressed as a linear coinbinatior 

1 ti. 

2 

of these past values of a^ and so its variance is proportional to o . Thus, 
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2 2 

Oj = a [1 + V(T)] . When dealiny with cumulative forecasts, the same aryunient 

holds. The error now is T ~ ^t T " ^t T show it has variance 

of the form o^[T + V^(T)] . 

2 

In practice, a is unknown. It may be estimated from the data, however, 
2 2 

using the fact that a [1 + V(l)] = , the variance of the one step ahead 

forecast error. Thus, we can estimate directly from the forecast errors, 

9 

and knowledge of V(l) yields an estimate of a . 

We derive V(T) and V^(T) from the variance-covariance matrix of the 

estimates of the individual components of the forecast i.e. m^, b^ and 

S. . , , k = 1 ,2 , . . . ,n . These variances are of some interest in their own 

t+k-n 

right for the construction of confidence intervals for the components, and we 
display them below. Details of the derivations are given in Appendix 1. 

= Var(m^.) = a^{a^[nd^ ^ ^ I (n-i)d.] + 2( 1-a^) (d^-d^) }o^ 

i = l 

= Var(S^^|^_^) = 2 q 2 ( l-OQ)^(dQ-d^)a^, k = l,2,...,n 
°mb ^ Cov(m^,b^) = aQaj^(2-Oj^) (dQ-d|^)o 



0|^k ^°''(^t’^t+k-n^ ~ ^^n-k'^^k ^k-1 " ^n-k+P® ’ 1,2,. ..,n . 

°mk "" ^°'^^"\’^t+k-n^ ^ ^'^n-k”^k^ ^^k-l~‘^n-k+l ^ ^ "" l,2,...,n 



o . • = a • • 
ij Ji 



(2dj_^ - - dj_.^j)o , 111 < j < n 



Before discussing the derivation of the values {d^ ,d^ , . . . ,d^ } we note the forms 
of V(T) and V^(T) . 
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Now, V(T)o‘^ = var(m^ + 



t t+k-n^ 



2 2 2 2 

“ 0 _n‘''Ta. +0. + 2 (oi + T0.+T0..) 
m D s mk 



mb 



bk^ 



Thus , 

V(T) = 2a^a^(dQ-d^)T^ 

+ 2aQa^[aQ(2-a^)(d|^-d^) + ^ 2 ^ ^'“o^ ^^n-k ' ^k ‘^k-l ' ^n-k+1^^''' 

+ 2aQa2(l-ocy)C(l-a^)(d^_^-d^) + ( d^_ ] 

+ 2oq( 1-a^ ) (dQ-d^) + 2a^( l-Oy)^(dQ-d^) + a^a^[ndy ^ ^ 1 (^"'')d^] (4) 

In the same way, using the fact that 

't,T' H*fT(TH)bj + j , 

k=l 

we can show that when T = pn + q, p ^ 0, 1 £ q _< n; 

V^(T) = Ll^ol + T^d+Do^^ + I T^d+D^o^l/o^ 

+ 2a^(l-aQ)2[(p+l)^(dQ-dq) + P^(dQ-d^_q) + p(p+l) (dq-d^+d^_q-dQ) ] 

+ aQa^a2(l-aQ)[(2p+l)(dy-d^) + ( d^_^-d^ ) ]T( T+1 ) 

+ 2aQCX2(l-c^Q)[p(2-ap(dQ-dJ + ( dQ-d^-d^+d^_^ ) 

^ -1 X ^'k-'n-k)^T (5) 

k = l 

Now the numbers dQ,d^,...,d^ are necessary to evaluate the variances. 

These (n+1) values are derived by means of an algorithm which is given and 
discussed in Appendix 2. It is a fairly straightforward algorithm due to 
Wilson (1979). It tias the property that it also checks the stability of the 
forecasting system. As we shall see in Section 3 this is most useful for tins 
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forecasting system. Before discussing any of the technical aspects of the 
procedure described here, we present some particular cases which have closed 
form solutions. 
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2. SPECIAL CASES 



2.1. Trend-free additive seasonal model. 

In this case the model assumed is a constant level plus addi- 
tive seasonals. The forecasting equations are obtained by removing the second 
equation from the original three (2), and dropping b^ ^ from the first. 

Here the values dQ,d^,...,d^ can be derived algebraically. The variance- 
covariance matrix of the components is now obtained from 



where 
Thus , 



h 



= o^/[(2- 



2 - n 2u 

^m - ’ 

“ oi2(l-aQ)[(n-l)o!Q + 

®mk " aQa2(l-Oy)h , 

‘^ij " ®ji " ■ °mk ’ 

aQ-a2+ag02) (00(0+02-0(002) ] 



02(1-00) ]h , 

1 £ k £ n , 

1 < i < j < n 



1 - (l-Of,)(l-Op) 

' 1 + (l-a^Hl-Cj) • 

It is interesting to note that for this model V(T) is independent of T . 

For the cumulative forecast error, we find that for T = pn + q, (p £ 0, 

1 £ q £ n). V^(T) = PT^ + o^( 1-O0)^( pT + pq + q) 

2 

where P = OgCnog + 02 ( 1 - 00 ) 3h/cr^ 



2.2. Non-seasonal Holt-Winters' system. 

The appropriate equations are obtained by removing the third equation 

from (2) and S. from the first. This is a well-known non-seasonal fore- 
t-n 

casting system and often used for linear trend models. It is a direct compet- 
itor of Brown's second order, or double, exponential smoothing which we will 
consider shortly. The underlying model is assumed to be a simple linear trend 

i.e. s is removed from (1). 
k 



7 



We obtain 



°mb = «o“i(2-ai)h^ , 

2 

where = a / (4-20^-01^0^ ) . 

We can deduce in the usual way that 

V(T) = [2ogO^T^ + 2oyOj^(2-Oj^)T + ( 2oQ+2oj^-3oQaj^ ) ]hj^/a^ (7) 

and 

V^(T) = [oqO^T^ + 4oQa^T + (4 o|j+4oj^-2oqOj^-OqO^) ]T^hj^/2o^ 

2.3. Brown's Double Exponential Smoothing. 

It is well-known that this system is equivalent to the Holt-Winters' 

2 

non-seasonal form above if we choose Oq = 1-3' and oj^ = (l-3)/(l+3) where 

3 = 1-0 and o is the smoothing constant. Thus, the corresponding results 

can be obtained by making these substitutions in the equations of Section 2.2. 

2 2 

above. The results for o^, and V(T) may be found in Brown (19v62) 

and Vj-(T) is given in Bowerman and O'Connell (1979). 

2.4. Continuously re-normal i zed seasonal factors. 

In the general additive seasonal model given by (1) it is usually 
n 

assumed that 'I s. = 0 . This is done to ensure some measure of independence 
k=l 

between the level of the data and the seasonal pattern. It is always recom- 
mended that when the correspond! ng forecasting system, as given by (2), is 
implemented, the seasonal factors sum to zero initially. What is to be done 
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thereafter is not so clear. Because of the revision equation for the seasonal 
factors they will no longer sum to zero after the first observation. We can 
renormalize the seasonal factors at any time by subtracting from each the 
average of the set, i.e. the most recent n . Whether we should do so or not 
is not clear. Some authors seem to recommend against it, e.g. Bowerman and 
O'Connell (1979); some regard it as an optional modification, e.g. Montgomery 
and Johnson (1976); some recommend renormalization once per season, e.g. 
Chatfield (1978); and some recommend continuous renormalization, i.e. after 
every revision of seasonal factors, e.g. Thomopoulos (1980). 

Our only interest here is in how such a procedure affects the forecast 
error variance. In the cases of occasional or purely seasonal renormalization 
the situation is very complex and we have nothing to say, except that the 
effect on the error variance appears to be small. In the case of continuous 
renormalization, however, the following result may be applied. An outline of 
the proof is given in Appendix 1. Continuous renormalization of the general 
forecasting system given by (2) yields exactly the same forecasts (and so 
forecast errors) as running the system without renormalization but replacing 
aQ,Oj,a 2 by ag»ct^ respectively, where a* = •> a 2 (l-o<g)/n , 

“l "" “o“l^“o ’ “2 "" ^ consequence, the appropriate 

covariance matrix elements and values of V(T) and Vj.(T) can be obtained 
as in Section 1 above by replacing Og, 0^,02 by ag,a^,a 2 respectively. The 
same holds true for the trend-free version discussed in Section 2.1. 
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3. TECHNICAL CONSIDERATIONS 



3.1. Infinite Error Variance and Stability. 

In smoothing systems in general the variance of the forecast error 
increases as the values of the smoothing constants increase. For example, in 
the case of the simple exponentially weighted moving average (SEWMA), 

Xt ” “X^ ^ , the variance of the forecast error is proportional to 

2/(2-a) which clearly increases with a in (0,1). The need for higher values 
of a in practice reflects the fact that the underlying level is changing 
rapidly. Consequently, a more responsive forecast is needed. Moreover, the 
inherent instability of the underlying model is reflected in the increase in 
forecast error variance. This variance is finite while a remains in the 
stability region of the system, i.e. (0,2) for the SEWMA. Note that the set 
of values from which a is usually chosen is a subset of the stability region 
of the system. Thus, a stable system always results for the SEWMA. This is 
also true for all the General Exponential Smoothing models. However, it is 
not true for the seasonal system under consideration here. 

It is a somewhat surprising and problematic fact that there are choices 
of the smoothing constants a^, , a 2 lying in the usual range, (0,1), which 

yield an infinite variance for the forecast error. The algorithm given in 
Appendix 2 checks for this possibility which indicates that the forecasting 
system is unstable. If such a situation arises clearly no meaningful confi- 
dence intervals can be constructed. More importantly, however, we must decide 
how to interpret this knowledge of the system's instability. The concept of 
stability for a system of difference equations is an important one but rarely 
discussed in the context of forecasting systems. Two useful exceptions to 
this are the papers by McClain (1974) and Brenner et al , (1968). 
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In essence, a forecasting system such as (2) above is stable if the 
influence of earlier observations decreases with the passage of time. Thus, 
the forecasts, or (equivalently) the forecast errors, are influenced more by 
recent observations than by those in the past. It is worth noting that this 
is also the essential philosophy of exponential smoothing systems, and indeed, 
most forecasting systems. 

The converse of this is that in an unstable system past observations have 

a constant or even growing influence on future forecasts. As an illustration, 

consider the SEWMA again. The forecast can be written as a weighted 

average of all past observations. At time t , the weight given to X is 

t “ k 

a(l-a) . Clearly, if a lies outside the stability region (0,2) this weight 
increases with k so that data in the most remote part have greatest influ- 
ence upon the forecast. Equally, if a = 2, all observations, however distant 
in time, make the same contribution to the forecast. 

The inescapable conclusion is that it would be extremely unwise to select 
smoothing constants which do not lie in the stability region of the system. 

On the one hand, we can tolerate a high (but finite) forecast error variance 
because this represents a trade-off between accuracy and robustness. In an 
effort to predict a model whose parameters are changing rapidly in time we 
may require a more responsive system. The cost of this is a correspondingly 
higher error variance. On the other hand, we can not tolerate an infinite 
error variance because this indicates an unstable forecasting system. Such a 
system violates exactly those assumptions which we hold most important to the 
generation of our forecasts. 

3.2. Practical consequences. 

It is obvious from the forms of V(T) in the two special cases (6) and 
(7) that the usual range of values of smoothing constants, i.e. (U,l), lies 
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within the stability regions for these systems. The problem arises for only 
the general seasonal model. The actual stability region is difficult to es- 
tablish in general here for it depends not only on Oq, and also 

on n , the length of the season. However, Gardner (1984) reports that Sweet 
has demonstrated numerically that the usual range for the three smoothing 
constants lies within the stability region for n up to four. For seasons 
longer than four periods it is no longer true; a result of importance for 
weekly and monthly data. 

The values of the smoothing constants are established in one of two ways 
in general. They may be selected intuitively by appealing to the ideas of 
required speed of response and constancy of the underlying model. Alterna- 
tively, they may be chosen as giving the best fit with respect to some cri- 
terion such as least squares. However the values are chosen, we should be 
very concerned if they do not lie in the systems 's stability region. 

If we have chosen them for rapid response we may allow their use for a 
brief period, perhaps when we initialize the forecasting system. As noted 
above, however, if we do not replace them with values in the stability region 
then the longer we forecast the more influential become the very first 
observations. If some best fitting criterion selects 'unstable' smoothing 
constants, i.e. ones leading to an infinite error variance, then a very real 
possibility is that the model is wrong, at least for part of the data. In 
particular, it may well be that a multiplicative seasonal model is called 
for. 

We may decide to retain the model but use smoothing constants in the 
stability region. These will be obtained by decreasing one or more of 
tty. Op . In general, if a^, ap 02 do not lie in the stability region 
it will be because one or more of them is too large. 
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4. SUMMARY 



The purpose of this paper has been to give a procedure for deriving the 
variance of the forecast error for Winters' Additive Seasonal System. In 
summary that procedure is as follows: 

(i) If the model has no trend or is non-seasonal the results are given 
in Section 2.1 and 2.2 respectively in closed form. 

(ii) From the chosen smoothing constants Oq, and the 

coefficients , W^,..., are derived as in Appendix 1. 

(iii) The W^. values are used to start the algorithm of Appendix 2 and 
yield d^j ,dj^ , . . . ,d^ . The algorithm simultaneously checks the 
stability of the forecasting system for the chosen values 
Oq, Op 02 . If n £ 4 the system will be stable. 

(iv) The values of d^. are substituted into the appropriate expressions 
for V(T) and V^(T) . 

Note that at step (iv) we are able to evaluate the variance of the errors for 
T-step ahead forecasts and T-step ahead cumulative forecasts. Confidence 
intervals for these future values can be derived from these variances. 

A discussion of how to calculate the correspond! ng results when the 
seasonal factors are renormalized after each observation is given. We have 
also discussed the (real) possibility that the error variance may be infinite 
and indicated its relationship with the stability of the forecasting system. 
Some recommendations are made about the interpretation of system instability. 
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APPENDIX 1 



The results of this paper are obtained usiny the fact that the forecast- 

iny system yiven by (2) has an equivalent ARIMA process. The nature of this 

equivalence is discussed in McKenzie (1984) and the equivalent ARIMA for (2) 

is yiven in McKenzie (1976). It has the form ( 1-B) ( 1 -b'^)X. = W(B)e. , where 

k k 

B is the backshift operator defined by B X = X. . and W(B) = 1 - ^ W.B 

k=l ^ 

where = 1 - = 1 - - 02 ( 1 - 0 ^) , = - (l-aQ)(l-a 2 ) ’ 

and W|^ = - for k = 2,3,...,n-l . 

Suppose now that W(B)6(B) = 1 where 6(B) = 'l 6,B^ . Usiny 

k=0 

equations (2), it can be shown that, for the purposes of the variance calcula- 
tion, the components of of the forecast can be expressed as infinite moviny 
averayes in {a^} . Thus: 



■"t ^ *-^i ■ ^^"“i^^i-i^*^ ^^t ’ 

00 , 

St = a2(l"Oo)Cl I ( » 

mi n( k ,n-l ) 

where Ci. = '[ <5. . . Further, the correspondi ny representation of S. . 

K i "“Q U’^K*"D 

n _ k 

can be obtained by writiny it as B , (k = l,2,.,.,n) . 

Since {a^} are independent random variables tlie variances and covariances 
can be obtained directly from these moviny averaye representations. Defininy 

00 

^ 1 / " 1 ~ yields the expressions yiven in the paper. 

K 1 1 K 

The alyebra is tedious and hardly illuminatiny and so is omitted. It can be 
obtained upon request from the author. 
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As regards the renormalization procedure, recall the equivalent ARIMA for 
(2) given above. If we introduce continuous renormal i zation as discussed we 
find there is still an equivalent ARIMA and it has the same form. Now, 
however, W(B) has Uq, replaced by o^, a^, . This ARIMA is equiv- 

alent to the system (2) with the starred smoothing constants. Hence, the 
result. Initially, the result is surprising. However, note that we are 
really dealing with two different decompositions of the seasonal factor. One 
form is appears in (2) and the other is where 

S*. is normalized to sum to zero and m* is the correspondi ngly adjusted 
t+k-n t 

level . 
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APPENDIX 2 



To obtain the values {d|^} note that, by definition, they are the variance 
(djj) and the first n autocovariances (d^ ,d 2 , . . . ,d^) of the autoreyressive 
process of order (n+1) given by W(B)Z^ = , where are independent 

random variables of zero mean and unit variance. Thus, the sequence {d|^} may 
be obtained by solving a suitable set of the Yule-Walker equations for this 
process. This procedure is discussed in McLeod (1975, 1977). 

From our point of view, however, a much superior approach is presented by 
Wilson (1979). The stability of the forecasting system corresponds to the 
stationarity of the autoregressive process {Z^} and can be tested routinely 
within the procedure. The algorithm is as follows: 

(i) define = Wj^ , k = 1,2,..., n+1; and t^^^ = 1 . 

(ii) apply the following equations in the given order for 
k = n+1 ,n,. . . ,2 

“k ' 1 - "Ik 

£ 0 , the system is unstable; stop now, 

^-l,i = ("^k,i ^ \,k ^^k,k-i)/\’ = l,2,...,k-l) 

^k-1 " ^k^'^k 

(iii) Dj = 1 - W^^^ 

if £ 0 , the system is unstable; stop now. 

d^ = t^/Dj^, and the system is stable, and the error variance finite, 

(iv) dj = d„ 



k-1 



'^k-l,i ^k-i ^\,k ^k 



, k — 2,3,... ,n 



Note that if n £ 4 the system is stable for the usual choice of parameters 
i.e. c (0,1), i = 0,1,2, and so we need not test Dj^ . 
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